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QUALITY  CONTROL  OF  METEOROLOGICAL  OBSERVATIONS 
Toke  Jayachandran  and  Richard  Franke 

Introduct ion  :  A  major  input  component  to  weather  prediction 
systems  is  the  observational  data,  typically  geopotential  heights 
and  winds  from  weather  observing  stations  around  the  world.  The 
Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS) 
produces  weather  forecasts  using  a  spectral  forecast  model. 
Preceding  the  forecast  is  a  multivariate  objective  analysis, 
followed  by  a  nonlinear  normal  mode  initialization  scheme.  The 
multivariate  objective  analysis  scheme  uses  the  weather 
observations  to  produce  analyzed  fields  with  minimum  statistical 
error.  As  is  to  be  expected,  some  of  the  observations  will  turn 
out  to  be  erroneous  for  various  reasons  such  as  the 
mis-calibration  of  the  measurement  devices,  failure  to  perform 
prescribed  data  corrections,  communications  errors  and  simple 
recording  errors  .  We  present  in  this  paper  a  quality  control 
methodology  to  screen  the  data  to  identify  those  observations  that 
appear  to  be  erroneous  in  a  statistical  sense.  It  is  then  up  to  a 
meteorological  analyst  or  a  decision  making  algorithm  to  examine 
the  tagged  observations  and  decide  on  a  course  of  action  -  delete 
the  miscreant  observations,  modify  the  data  if  appropriate  or 
include  the  observations,  as  is,  for  NOGAPS  analysis.  The 
procedures  can  also  be  used  to  monitor  the  reliability  of  the 
observing   stations  and  to  alert  them  of  any  mispract ices . 

Data  Quality  Tests:  The  methodology  proposed  here  is  an 
adaptation  of  a  statistical  test  for  outliers  (Dixon,  1950  [1]). 
The  first  step  in  the  application  of  this  test  is  to  partition  the 
monitoring  stations  into  small  "contiguous  groups" .  The  recorded 
observations  within  any  single  group  are  expected  to  be  "similar" 
because  the  weather  patterns  at  these  monitoring  sites  are 
comparable  due  to  the  proximity  of  the  stations  to  each  other. 


One  approach  towards  identifying  such  a  group  of  stations  is  the 
following.  Select  the  latitude  and  longitude  for  a  target 
location  on  the  surface  of  the  earth  (this  could  be  the 
coordinates  of  a  reporting  station  of  interest)  and  choose  an 
approximately  square  region  by  specifying  a  range  for  the 
latitude,  say  5°  above  and  below  the  target,  and  within  a  similar 
distance  in  the  longitudinal  direction.  The  observational 
increments  (differences  between  the  reported  measurements  and  the 
forecast  values)  for  all  of  the  stations  within  the  designated 
area  are  then  subjected  to  a  sequence  of  outlier  tests. 

Assume  that  the  number  of  reporting  stations  within  the  target 
area  is  N  and  the  data  from  the  most  recent  k  days  (say  30  days) 
is  to  be  examined  for  quality.   Let   Xi!  ,  Xi2  ,  .  .  .  ,  Xik    be 

the  observation  increments  for  the  ith  station,  i  =  1,  2,  .  .  .  , 
N;  let  Mx  ,  M2  ,  .  .  .  ,  MN  be  the  means  of  the  increments 
(averaged  over  the  k  days  for  each  of  the  N  stations)  and  Si,  S2  , 
.  ,  SN  the  standard  deviations  (over  the  k  days)  for  the  N 
stations . 

The  first  test  we  propose  is  to  compare  the  daily  observational 
increments   Xij  ,  X2j  ,  .  .  .  ,  XNj   for  the  jth   day,  and  check  for 

"outliers".    Let   Yx  <  Y2  <  .  .  .  <  YN  be  the  data  rearranged  in 

increasing  order.   Compute 


rn  =  (YN  -  YN_!)  /  (YN  -  Y2)  ; 

if  r:i  exceeds  the  tabulated  value  in  Table  A  (extracted  from 
Dixon, 1950) ,  corresponding  to  a  chosen  statistical  significance 
level  ft  (say  .05),  conclude  that  the  observed  datum  for  the 
station  corresponding  to  YN  is  an  outlier.   This  implies  that  the 


observation  increment   YN  is  too  large  relative  to  the  increments 

for  the  other  stations  in  the  region.  That  is,  the  observed 
geopotential  height  appears  to  be  too  high  (large)  when  compared 
to  the  forecast .   Now  compute 


r'u  =   (Y2  -Yi)  /  (YN_!  -  Yi) 

and  conclude  that  Y:  is  an  outlier  if  r'u    exceeds  the  tabulated 

value  in  Table  A;  if  this  happens,  the  implication  is  that  the 
observed  geopotential  height  is  too  low  (small)  compared  to  the 
forecast.  These  tests  can  be  performed  for  each  of  the  days  for 
which  data  is  available.  The  meteorological  analyst  or  the 
decision  making  algorithm  will  then  have  to  choose  a  course  of 
action,  based  on  some  established  decision  rule,  for  each  of  the 
observations  failing  the  test.  The  two  tests  were  applied  to 
actual  data  from  three  different  areas  of  the  globe  viz.,  US, 
Europe,  and  China.  In  each  area  we  identified  a  location  by 
specifying  its  latitude  and  longitude  and  then  selected  all  the 
stations  within  4°  of  this  location.  The  results  of  the  daily 
tests  are  presented  in  Tables  la,b,c.  The  I.D.  numbers  of  the 
stations  included  for  the  test  are  listed  in  the  first  column. 
The  entries  (there  are  two  entries  in  each  of  the  16  columns)  in 
the  table  represent  the  number  of  days  in  October  *  90  on  which  the 
stations  failed  the  outlier  test  at  the  low  end  (rn'  test)  and  the 
test  at  the  high  end  (rn  test)  using  the  data  (geopotential  height 

differences)  for  16  pressure  levels  (1000,  850,  700,  500,  400, 
300,  250,  200,  150,  100,  70,  50,  30,  20,  15,  10  mb) .   A  dot  in  any 
position  implies  no  failures;  an  asterisk  indicates  insufficient 
data  to  perform  the  tests.    Such  a  table  can  be  used  in  making 
quality  judgments  on  the  performance  of  the  monitoring  stations. 

Another  test  that  can  also  be  used  to  monitor  the  performance 


quality  of  weather  stations  is  to  perform  the  two  outlier  tests  on 
the  averages  Mi  ,  M2  ,  .  .  .  ,  MN  .  The  data  for  the  tests  are 
these  averages  rearranged  in  increasing  order.  Tables  2a,b,c 
present  the  same  type  of  information  as  Tables  la,b,c  except  that 
this  time  the  entries  are  a  1  or  a  dot;  a  1  signifies  that  the 
average  for  a  station  failed  the  outlier  test  while  a  dot  would 
signify  that  the  station  passed  the  test. 

A  third  outlier  test,  that  can  be  viewed  as  a  test  for 
consistency,  is  to  check  for  an  outlier  in  the  station  standard 
deviations  Si  ,  S2  ,  .  .  .  ,  SN  .  If  Yi  <  Y2  <  .  .  .  <  YN  are 
these  standard  deviations  in  ordered  form,  compute 

r10  =  (YN  -  YN_!)  /  (YN  -  Yx) 

and  conclude  that  the  standard  deviation  corresponding  to  YN  is 

statistically  too  large  relative  to  the  others,  if  ri0  exceeds  the 

critical  value  in  Table  B  (extracted  from  Dixon,  1950)  at  a 
significance  level  ii  (say  .05)  .  Tables  3a,b,c  are  similar  to 
Tables  1  and  2  using  the  station  standard  deviations,  with  the 
difference  that  only  the  high  end  outlier  test  is  used. 

We  have  written  a  computer  program,  using  Fortran,  that  will 
select  (1)  the  reporting  stations  within  a  specified 
(approximately  square)  region  around  a  specified  location  and  (2) 
the  appropriate  data  for  these  stations  from  the  Naval 
Oceanographic  and  Atmospheric  Research  Laboratory  (NOARL)  files 
[2],  and  automatically  performs  all  of  the  tests  described  above. 
Tables  1,  2  and  3,  described  above,  are  the  actual  outputs  of  this 
program.  The  program  also  generates  additional  information  such 
as  the  number  of  times  the  daily  outlier  tests  were  performed 
(Table  4)  and  the  number  of  days  for  which  data  was  available 
(Table  5),  during  the  period  specified  for  the  test.   The  numbers 


in  Table  4  will  sometimes  be  smaller  than  the  corresponding 
entries  in  Table  5;  this  happens  whenever  the  number  of  data 
points  available  for  the  test  is  smaller  than  the  minimum  sample 
size  of  4.  Table  6  is  an  example  of  another  type  of  output  of  the 
program;  for  a  selected  station  the  table  provides  the  daily 
test  history,  over  the  16  pressure  levels.  A  "PASS"  indicates 
that  the  station  passed  all  the  daily  tests;  a  "LFAIL"  signifies 
that  the  observational  increment  is  too  low  as  compared  to  the 
numbers  for  the  other  stations  in  the  region,  on  at  least  one  day 
of  the  testing  period.  A  "FAILH"  denotes  that  the  increment  for 
at  least  one  day  was  deemed  to  be  high.  Similar  tables  displaying 
the  results  of  the  means  and  standard  deviations  tests  are  also 
produced.  In  the  "TOTAL"  column  a  "PASS"  means  that  the  station 
passed  all  the  tests  on  all  the  days  and  a  "FAIL"  implies  the 
converse. 
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STATION  /  Lvl  1 


9  10  11   12   13  14   15  16 


72456 
72357 
72349 

2  . 

..2.5.2.1 
.11 

.1.1 
1  .  2 

.113.  .  .2.  .  .2.  .112** 
.1.1.  1.2.  ..2.  1.1.** 
1 *  * 

72340 

11 1  .  *  * 

72260 

.2 

.  1  . 

121 1.1.  ..2** 

72247 
72239 

1.11.1... 
1  .  .  . 

.  .  1 

...1...  1.1.1.3.5.** 

,  ,  ,**************** 

72257 

************* 

******************* 

Table  la:   Failures  in  daily  outlier  test  (.  =0,  *  =  no  test) 

U.S.  stations  near  (35  ,-95  ) 

10/01/90  to  10/31/90  at  00Z 


STATION  /  Lvl  1   2   3   4   5   6   7   8   9   10  11   12   13  14   15  16 

26702  .1.1.1.3.1  1.1.1.1.2.1.. 

26406  1  .  1  .  2  .  1  .  1  .  1  .  .  3  .  .  .  1 1  . 

12425  .  .  1  .  .  1 1  .  .  1 1.21211.... 

12374  1.1 1  .  1  .  1  .  1  .  .  .  2  .  2  .  1  .  . 

12330  ..11 1  ....  - 1 1.1**** 

12120  1 **** 

10437  .1.1 1 ********** 

10338  

10046  1.1.2  ********** 

10035  1  .  .  . 

9393  1 1  .  .  . 

9184  1  ...  1 

6181  

2591  1.2.1.1 **** 

2527  1.1 **** 


Table  lb:   Failures  in  daily  outlier  test  (.  =0,  *  =  no  test) 

European  stations  near  (55  ,15  ) 

10/01/90  to  10/31/90  at  00Z 


STATION  /  Lvl  1   2   3   4   5   6   7   8   9   10   11   12   13   14   15  16 

57679  1 1.1.2 ** 

57516  **.l.l..l  1. .12.1.1.1.1** 

57494  .11.2.11. .1.1.1  ** 

57461  ..**..l 1 ** 

57447  .3** ** 

57328  ****,  .1.1  ****************** 

57178  .  .  *  *  1  .  .  2.  11.  .112.  3.  1.1 ** 

57127  ..21.  .11.1 1 1** 

57083  ...1...1  1  1...** 

57036  .41.111.  .11 1 ** 

57245  ****  ********** 

57633  ******   i    ,    .  ,********************** 

Table  lc:   Failures  in  daily  outlier  test  (.  =0,  *  =  no  test) 

Chinese  stations  near  (32  ,110  ) 

10/01/90  to  10/31/90  at  12Z 
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STATION  /  Lvl  1 


9  10  11   12   13  14  15  16 


72456 
72357 
72349 
72340 
72260 
72247 
72239 
72257 


....1.1 

..........................  i  !  i  .  .  '. 

******************************** 


Table  2a:   Failures  in  the  means  test 


=0,  *  =  no  test) 

-O 

) 
10/01/90  to  10/31/90  at  00Z 


U.S.  stations  near  (35  ,-95  ) 


STATION  /  Lvl  1 


9  10  11   12  13  14  15  16 


26702 

26406 

12425 

12374 

12330 

12120 

10437 

10338 

10046 

10035 

9393 

9184 

6181 

2591 

2527 


*  * 

,  ,  .  ,****** 

********** 
*  *  *  * 


Table  2b:   Failures  in  the  means  test  (.  =0,  *  =  no  test) 

European  stations  near  (55  ,15  ) 

10/01/90  to  10/31/90  at  00Z 


STATION  /  Lvl  1   2   3   4   5   6   7   8   9  10  11   12   13  14   15  16 

57679      *  * 

57516      **.l 1 1.1.1.. 

57494      

57461      ..** ** 

57447      .  .  *  *  

57328      ****  i    t  #  .****************** 

57178      ..** ** 

57127      

57083      *  * 

57036      *  * 

57245      ****  ********** 

57633  ******************************** 


Table  2c:   Failures  in  the  means  test  (.  =0,  *  =  no  test 

Chinese  stations  near  (32  ,110  ) 

10/01/90  to  10/31/90  at  00Z 


8 


STATION  /  Lvl  1 


9   10   11   12   13   14   15   16 


72456 
72357 
72349 
72340 
72260 
72247 
72239 
72257 


1   1 

*  *   * 

*  *   * 


Table  3a:   Failures  in  the  standard  deviations  test 

U.S.  stations  near  (35  ,-95  ) 
10/01/90  to  10/31/90  at  00Z 


=  0,  *  =  no  test! 


STATION  /  Lvl  1 


8   9  10  11   12   13  14   15   16 


26702 

26406 

12425 

12374 

12330 

12120 

10437 

10338 

10046 

10035 

9393 

9184 

6181 

2591 

2527 


1 

* 

1 ** 

***** 

*   * 

*   * 


Table  3b:   Failures  in  the  standard  deviations  test  (.  =0,  *  =  no  test) 

European  stations  near  (55  ,15  ) 
10/01/90  to  10/31/90  at  00Z 


STATION  /  Lvl  1    2   3   4   5   6   7   8   9   10   11   12   13   14   15   16 

57679  * 

57516  *   

57494  

57461  .    * * 

57447  .    *   

57328  ** ********* 

57178  .   * * 

57127  

57083  * 

57036  * 

57245  ** ***** 

57633  **************** 


Table  3c:   Failures  in  the  standard  deviations  test  (.  =0,  *  =  no  test) 

Chinese  stations  near  (32  ,110  ) 
10/01/90  to  10/31/90  at  00Z 


STATION  /  Lvl  1 


9  10  11   12  13  14   15  16 


72456 

31 

30 

31 

31 

31 

31 

31 

31 

31 

31 

30 

29 

27 

24 

23 

0 

72357 

31 

30 

31 

31 

31 

30 

31 

31 

31 

30 

30 

30 

29 

26 

22 

0 

72349 

31 

30 

31 

31 

30 

31 

30 

30 

30 

30 

30 

29 

28 

28 

22 

0 

72340 

31 

29 

31 

31 

30 

30 

30 

30 

30 

30 

29 

28 

28 

26 

22 

0 

72260 

30 

29 

30 

30 

30 

30 

30 

29 

29 

29 

29 

29 

28 

26 

23 

0 

72247 

31 

28 

31 

31 

31 

31 

31 

31 

31 

31 

29 

29 

29 

27 

23 

0 

72239 

1 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

72257 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Table  4:   Number  of  days  daily  outlier  test  was  performed 

U.S.  stations  near  (35  ,-95  ) 

10/01/90  to  10/31/90  at  00Z 


STATION  /  Lvl  1 


9  10  11   12   13  14   15  16 


72456 

31 

30 

31 

31 

31 

31 

31 

31 

31 

31 

30 

29 

28 

25 

24 

11 

72357 

31 

30 

31 

31 

31 

30 

31 

31 

31 

30 

30 

30 

30 

28 

25 

4 

72349 

31 

30 

31 

31 

30 

31 

30 

30 

30 

30 

30 

29 

28 

28 

23 

6 

72340 

31 

29 

31 

31 

30 

30 

30 

30 

30 

30 

29 

28 

28 

27 

24 

2 

72260 

30 

29 

30 

30 

30 

30 

30 

29 

29 

29 

29 

29 

29 

28 

27 

3 

72247 

31 

28 

31 

31 

31 

31 

31 

31 

31 

31 

29 

29 

29 

27 

23 

6 

72239 

1 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

72257 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Table  5 

:   Numbe 

r  of 

dai 

ly  r 

epor 

ts 

U.S.  stations  near  (35  ,-95  ) 
in/m  /op  4-^  10/91  /on  o+  nn7 


STATION 

72456: 

LAT , LON 

39, 

-95 

LEVEL  #DAYS 

DAILY 

TESTS 

MEANS 

TESTS 

STD  TESTS 

TOTAL 

1 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

2 

30 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

3 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

4 

31 

PASS 

*FAILH 

PASS 

PASS 

PASS 

*FAIL* 

5 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

6 

31 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

*FAIL* 

7 

31 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

*FAIL* 

8 

31 

LFAIL* 

*FAILH 

PASS 

PASS 

PASS 

*FAIL* 

9 

31 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

10 

31 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

11 

30 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

12 

29 

PASS 

PASS 

PASS 

PASS 

PASS 

PASS 

13 

28 

LFAIL* 

PASS 

PASS 

PASS 

PASS 

♦FAIL* 

14 

25 

PASS 

*FAILH 

PASS 

PASS 

PASS 

♦FAIL^ 

15 

24 

LFAIL* 

*FAILH 

PASS 

PASS 

PASS 

♦FAIL^ 

16 

11 

NOTEST 

NOTEST 

PASS 

PASS 

PASS 

PASS 

Table  6:   Summary  report  for  station  72456 
10/01/90  to  10/31/90  at  00Z 
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Appendix:   Program  listing 

C       THIS  PROGRAM  IS  FOR  CHECKING  THE  STANDARD  FORECAST  DEVIATION 

C       FILES  FOR  OUTLIERS  IN:   (1)   DAILY  READINGS,  (2)   MEANS,  AND 

C       (3)   STANDARD  DEVIATIONS  OF  STATIONS  WITHIN  A  TARGET  DISTANCE 

C       OF  A  GIVEN  LATITUDE  AND  LONGITUDE. 

C 

C       WRITTEN  BY  RICHARD  FRANKE 

C  DEPARTMENT  OF  MATHEMATICS 

C  NAVAL  POSTGRADUATE  SCHOOL 

C  MONTEREY,  CA  93943 

C  0083P@CC.NPS.NAVY.MIL 

C  408/646-2758 

C 

C       DATE  OF  BEGINNING:   6/21/91 

C       DATES  OF  CHANGES:    6/21/91 

C       MODIFIED  TO  USE  SYMBOLS  FOR  OUTPUT:   6/21/91 

C       MINOR  PRETTYING  UP:   6/24/91 

C       MODIFIED  TO  NOTE  NUMBER  OF  TESTS  AND  SUMMARY  BY  STATION:  6/27/911 

C 

PARAMETER  (MX=152 , NS=15 ) 

REAL  FE ( MX , NS , 16 ) , FER ( MX , NS ) , XMN ( NS ) , STD( NS ) , FED( NS ) 

I NTEGER  YRMD ( MX , NS ) , NSTAT ( NS ) , KNT ( NS ) , NSD( NS ) , 

1  NDLF(NS,16),NDHF(NS,16),NSDF(NS,16),NMNLF(NS,16),NMNHF(NS,16), 

2  NSS(NS),YMD(MX,NS),NDG(NS,16),NDDT(NS,16),LATS(NS),LONS(NS) 
CHARACTER* 15  IUNIT.OUNIT 

CHARACTER*6  FLAGS ( -1 : 1 ) , FLAGSL( -1 : 1 ) , FLAGSH( -1 : 1 ) 

CHARACTER* 1  SY(-1:36) 

CHARACTER*38  SYY 

EQUIVALENCE  ( SY , SYY ) 

DATA  SYY/ ' * . 123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ& ' / 

DATA  FLAGS/ ' NOTEST ' , '  PASS  ' , ' *FAI L* ' / , 

1  FLAGSL/' NOTEST' , '  PASS  ','LFAIL*'/, 

2  FLAGSH/ ' NOTEST ' , '  PASS  ' , ' *FAI LH ' / 
C 

C       GET  DATA  INPUT  FILE  NAME  AND  OUTPUT  FILE  NAME 
C 

WRITE(*,9) 
9   FORMAT ('  NAME  THE  INPUT  FILE') 
READ(*,10)IUNIT 
10   FORMAT (Al 5) 

OPEN(4,BLANK='ZERO' ,FILE=IUNIT,MODE='READ' ,STATUS= 'OLD' ) 
WRITE(*,21) 
21   FORMAT ('  NAME  THE  OUTPUT  FILE') 
READ(*,10)OUNIT 

OPEN(8,FILE=OUNIT,MODE='WRITE' ,STATUS= 'NEW' ) 
C 

C       GET  THE  INPUT 
C 

100   WRITE(*,1) 

1    FORMAT ('  INPUT  THE  TARGET  LATITUDE,  LONGITUDE  AND  THE  RANGE',/ 
1  '  -  FREE  FIELD  -  RANGE  =  0  STOPS  EXECUTION',/) 
READ(*,*)LAT,LON,LLR 
IF(LLR.LE.0)STOP 


J  i 


c 

C       STRETCH  OUT  IN  THE  LONGITUDINAL  DIRECTION 
C 

LTR=LLR/COS ( FLOAT ( LAT ) /5  7 . 3 ) 
C 

C       START  INPUT  UNIT  4  AT  THE  BEGINNING 
C 

REWIND  4 

MS  =  0 
C 

C       READ  A  HEADER 
C 

200   READ(4,2,END=300)NST,ILAT,ILON,ITM,NREC 

2  FORMAT ( BZ , 3X , 1 7 , 4X , 21 5 , 8X , 1 2 , 1 5 ) 
IF(ABS(ILAT-LAT).LE.LLR  .AND.  ABS( ILON-LON) . LE. LTR)  THEN 

C 

C         IT'S  IN  THE  RECTANGLE,  GET  THE  DATA 

C 

MS=MS  +  1 

NSTAT(MS)=NST 

LATS(MS)=ILAT 

LONS(MS)=ILON 

READ(4,3)((FE(I,MS,J),J=1,16),YRMD(I,MS),I=1,NREC) 

3  FORMAT(BZ,16F4.0,I8) 
ELSE 

C 

C         OR  SKIP  IT 

C 

READ( 4,4) (JUNK, 1=1, NREC) 
ENDIF 

4  FORMAT(Il) 
C 

C       GO  BACK  FOR  MORE  IF  WE  DON'T  HAVE  THE  MAX 
C 

IF(MS.LT.NS)  GO  TO  200 
PRINT  7 
7   FORMAT ('  STOPPING  BEFORE  END') 

300  CONTINUE 
IF(MS.LT.4)  THEN 

WRITE(*,5)MS 
C 

C         NOT  ENOUGH  DATA,  START  OVER 
C 

5  FORMAT ('  ONLY ',12,'  STATIONS  FOUND  -  START  OVER') 
GO  TO  100 

ELSE 

WRITE(*,6)MS,(NSTAT(N),N=1,MS) 

6  FORMAT (1 5,'  STATIONS  FOUND: ',/( 110) ) 
ENDIF 

C 

C       GET  THE  DATE  RANGE 

C 

301  WRITE(*,8) 
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8   FORMAT ('  INPUT  THE  BEGINNING  AND  ENDING  DAYS  (>=1,  <  =  152)V 

1  '  BEGINNING  =  ZERO  GOES  TO  NEW  LATITUDE/LONGITUDE'/ 

2  '  BEGINNING  =  NEGATIVE  ENDS  RUN') 
READ(*,*)IDB,IDE 

IF(IDB.LT.0)STOP 
IF(IDB.EQ.0)GO  TO  100 
C 

C       OK,  LET'S  DO  THE  JOB 
C 

C       ZERO  OUT  FLAG  FIELDS 
C 

DO  304  L=l,16 
DO  303  N=1,NS 
NDLF(N,L)  =  0 
NDHF(N,L)  =  0 
NSDF(N,L)  =  -1 
NMNLF(N,L)  =  -1 
NMNHF(N,L)  =  -1 
NDG(N,L)  =  0 
NDDT(N,L)  =  0 

303  CONTINUE 

304  CONTINUE 
C 

C       LOOP  OVER  ALL  LEVELS 
C 

DO  400  L=l , 16 
KS=0 
C 

C         KS  COUNTS  STATIONS  WITH  A  HISTORY  OF  2  OR  MORE  DAYS 
C 

DO  305  N=1,MS 
KNT(N)  =  0 
XMN(N)  =  0. 
C 

C  KNT  COUNTS  HOW  MANY  DAYS 

C 

305  CONTINUE 

DO  310  I=IDB,IDE 
ND=0 

DO  308  N=1,MS 
C 

C  ND  COUNTS  THE  NUMBER  OF  STATIONS  WITH  GOOD  DATA  TODAY 

C  FED  SAVES  IT  FOR  THE  TEST 

C  NSD  TELLS  WHICH  STATIONS  HAVE  GOOD  DATA 

C  FER  SAVES  IT  ALL  FOR  MEAN  AND  STD  CALCS 

C  YRMD  SAVES  THE  DATE  INFORMATION  (NOT  USED  PRESENTLY,  6/21 

C 

FED(N)=0. 

IF(FE(I,N,L).LT.999)  THEN 
NDG(N,L)=NDG(N,L)+1 
KNT(N)  =  KNT(N)  +  1 
ND=ND+1 

FED(ND)=FE(I,N,L) 
NSD(ND)=N 
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FER(KNT(N),N)    =   FE(I,N,L) 
YMD(KNT(N),N)    =   YRMD(I,N) 
XMN(N)=XMN(N)+FER(KNT(N),N) 
ENDIF 

308  CONTINUE 
C 

C  TEST  CURRENT  DAY'S  DATA 

C 

CALL  R11P05(FED,ND,NB,IB,NT,IT,IER) 
IF(IER.LE.O)  THEN 

IF(NB.GT.O)NDLF(NSD(IB),L)  =  NDLF(NSD( IB) ,L)  +  1 
IF(NT.GT.O)NDHF(NSD(IT),L)  =  NDHF(NSD( IT) , L)  +  1 
DO  309  N=1,ND 

NDDT(NSD(N),L)  =  NDDT(NSD(N) ,L)  +  1 

309  CONTINUE 
ENDIF 

310  CONTINUE 

DO  315  N=1,MS 

IF(NDDT(N,L).EQ.O)  THEN 
NDLF(N,L)  =  -1 
NDHF(N,L)  =  -1 
ENDIF 
315     CONTINUE 

DO  330  N=1,MS 

IF(KNT(N).GT.1)THEN 
C 

C  CALCULATE  MEAN  AND  ST.  DEV  IF  MORE  THAN  ONE  DAY  OF  DATA 

C 

KS=KS+1 
NSS(KS)  =  N 
XMN(KS)=XMN(N)/KNT(N) 
S=  0. 
DO  320  I=1,KNT(N) 

S=S+(FER(I,N)-XMN(KS))**2 
320         CONTINUE 

STD(KS)  =  SQRT(S/(KNT(N)-1)) 
ENDIF 
330     CONTINUE 
C 

C         RUN  TEST  ON  MEAN  VALUES  FOR  THE  STATIONS 
C 

CALL  R11P05(XMN,KS,NB,IB,NT,IT,IER) 
IF(IER.EQ.O)THEN 
DO  340  N=1,KS 

NMNLF(NSS(N),L)  =  0 
NMNHF(NSS(N),L)  =  0 
340  CONTINUE 

IF(NB.NE.O)  NMNLF(NSS(IB),L)=1 
IF(NT.NE.O)  NMNHF(NSS(IT),L)=1 
ENDIF 
C 

C         RUN  TEST  ON  STANDARD  DEVIATIONS  OF  THE  STATIONS 
C 

CALL  R10P05(STD,KS,NT,IT,IER) 
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IF(IER.EQ.O)THEN 
DO  360  N=1,KS 

NSDF(NSS(N),L)  =  0 
360       CONTINUE 

IF(NT.GT.O)NSDF(NSS(IT),L)  =  1 
ENDIF 
400   CONTINUE 
MDATB=0 
MDATE=0 
DO  410  N=1,MS 

I F ( MDATB . EQ . 0 ) MDATB= YRMD ( I DB , N ) 
I F ( MDATE . EQ . 0 ) MDATE= YRMD ( I DE , N ) 
410   CONTINUE 

I F ( MDATB . EQ . 0 ) MDATB=MDATE- ( I DE- I DB ) 
IF(MDATB.LT.0)MDATB=IDB 
I F ( MDATE . EQ . 0 )MDATE=MDATB+ ( I DE- I DB ) 
C 

C       OUTPUT  RESULTS  FOR  VARIOUS  TESTS 
C 

WRITE(8,ll)MDATB,MDATE,ITM,LAT,LON,LLR,LTR 

11  FORMATCl   RESULTS  OF  THE  OUTLIER  TESTS  FROM', 

1  18, '  TO', 18,'  AT',I3,'Z'// 

2  '    NUMBER  OF  FAILURES  BY  LEVEL  AND  STATION:  .=0,  *=NOTEST' 

3  //'    TARGET  LATITUDE/LONGITUDE  IS', 215,',  RANGE: ' , 12 , 'X' , 12 

4  //'    RESULTS  OF  THE  DAILY  OUTLIER  TEST  (LOW  FAIL/HIGH  FAIL)' 
WRITE(8,12)((L),L=1,16) 

12  FORMAT!/'  STATION  /  Lvl',1614/) 
WRITE(8,13)(NSTAT(N), 

1      (SY(NDLF(N,L)),SY(NDHF(N,L)),L=1,16),N=1,MS) 

13  FORMAT(I8,5X,32A2) 
WRITE(8,20)IDE-IDB+1 

20    FORMAT (//'    NUMBER  OF  DAYS  DAILY  OUTLIER  TEST  WAS  PERFORMED' 
1     '  OUT  OF ',14) 

WRITE(8,12)((L),L=1,16) 

WRITE(8,18)(NSTAT(N),(NDDT(N,L),L=1,16),N=1,MS) 
18   FORMAT(I8,4X,16I4) 

WRITE(8,14) 

14  FORMAT (//'    RESULTS  OF  THE  MEANS  OUTLIER  TEST  (LOW  FAIL/HIGH' 
1         '  FAIL)') 

WRITE(8,12)((L),L=1,16) 
WRITE(8,13)(NSTAT(N), 
1      (SY(NMNLF(N,L)),SY(NMNHF(N,L)),L=1,16),N=1,MS) 
WRITE(8,15) 

15  FORMAT(//'    RESULTS  OF  THE  STANDARD  DEVIATION  OUTLIER  TEST') 
WRITE(8,12)((L),L=1,16) 

WRI TE ( 8 , 1 6 ) ( NSTAT ( N ) , ( SY ( NSDF ( N , L ) ) , L= 1 , 1 6 ) , N= 1 , MS ) 

16  FORMAT(I8,4X,16A4) 
WRITE(8,17)IDE-IDB+1 

17  FORMAT (//'    NUMBER  OF  DAYS  REPORTING  OUT  OF ',14) 
WRITE(8,12)((L),L=1,16) 

WR I TE ( 8 , 1 8 ) ( NSTAT ( N ) , ( NDG ( N , L ) , L= 1 , 1 6 ) , N= 1 , MS ) 
C 

C       OUTPUT  SUMMARY  BY  STATION 
C 
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WRITE(8,22) 

22  FORMAT(//'  ************', 

2  //'  SUMMARY  BY  STATION  FOLLOWS'/) 
DO  500  N=1,MS 

WRITE(8,23)NSTAT(N),LATS(N),LONS(N) 

23  FORMAT (//'  


> 


2  /'  STATION', 16,':   LAT,LON  =  ',14, ',',14/ 

3  /'  LEVEL  #DAYS  DAILY  TESTS         MEANS  TESTS    ', 

4  '  STD  TESTS  TOTAL'/) 
DO  420  L=l,16 

C 

C  SET  CODES  FOR  TESTS:   -1  IS  NOTEST,  0  IS  PASS  ALL,  1  IS  FAIL 

C  ONE  OR  MORE 

C  FOR  THE  TOTAL  TEST,  NOTEST  MEANS  NO  TESTS  PERFORMED, 

C  PASS  MEANS  PASS  ALL  TESTS  PERFORMED  (IF  ANY) 

C  FAIL  MEANS  FAILED  AT  LEAST  ONE 


C 

NFL=MIN(NDLF(N,L),1) 
NFH=MIN(NDHF(N,L),1) 
NMFL=MIN(NMNLF(N,L),1) 
NMFH=MIN(NMNHF(N,L),1) 
NSF=MIN(NSDF(N,L),1) 
NTF=MAX ( NFL , NFH , NMFL , NMFH , NSF ) 
WRI TE ( 8 , 24 ) L , NDG ( N , L ) , FLAGSL ( NFL ) , FLAGSH ( NFH ) , 
1     FLAGSL ( NMFL ) , FLAGSH ( NMFH ) , FLAGS ( NSF ) , FLAGS ( NTF ) 
420     CONTINUE 

24     FORMAT(2I6,4X,A6,3X,A6,5X,A6,3X,A6,5X,A6,5X,A6) 
500   CONTINUE 

WRITE(8,19) 
19   FORMAT(//'  ******************', 

I'*********************************************************' ///// ) 
GO  TO  301 
END 

SUBROUTI NE  R10P05 ( X , N , NOUTT , I T , I ER ) 

PARAMETER  (MAX=30) 

DIMENSION  X(N),SX(MAX),ISX(MAX) 

DIMENSION  CV(30) 
C 

C     THIS  SUBROUTINE  DOES  THE  RIO  TEST  FOR  OUTLIERS  AT  THE  .05  CRITICAL 
C     TEST  LEVEL. 
C 

C     THE  CV  VALUES  ARE  THE  CRITICAL  VALUES  FROM  TABLE  V. 
C 

DATA  CV/2*0,. 941,. 765,. 642,. 560,. 507,. 468,. 437,. 412,. 392,.  376, 

1  .361,. 349,. 338,. 329,. 320,. 313,. 306,. 300,. 295,. 290,. 285,. 281, 

2  .277, .273, .269, .266, .263, .260/ 
C 

C  INPUT  VALUES  ARE:   X  -  THE  DATA  TO  BE  TESTED 

C  N  -  NUMBER  OF  DATA  POINTS 

C 

C  OUTPUT  VALUES  ARE:   NOUTT  -  FLAG  FOR  LARGEST  VALUE 


16 


c 

c 
c 
c 
c 
c 
c 
c 
c 
c 


0  IF  OK,  1  IF  NOT 
IB    -  LOCATION  IN  X  ARRAY 
IER   -  ERROR  RETURN  INDICATOR 

0  IF  ALL  OK 

1  IF  N  OUT  OF  RANGE  4  TO  MAX,  INCLUSIV 

NO  TEST  PERFORMED 


CHECK  FOR  SUITABLE  N 

IER  =  0 

NOUTT  =  0 

IF(N.LT.3.0R.N.GT.MAX)THEN 

IER  =  1 

RETURN 
ENDI  F 
DO  100  1=1, N 

SX(I)  =  X(I) 

ISX(I)  =  I 
100  CONTINUE 

SORT  INTO  INCREASING  ORDER 

CALL  SHSORT(SX,ISX,N) 

RT  =  0. 

IF(SX(N)-SX(1)  .NE.  0.)  RT  =  (SX(N)  -  SX(N-1 )  )/(SX(N)  -  SX(D) 

IF(RT.GT.CV(N))THEN 

NOUTT  =  1 

IT  =  ISX(N) 
ENDIF 
RETURN 
END 


SUBROUTI NE  Rl 1 P05 ( X , N , NOUTB , I B , NOUTT , I T , I ER ) 
DIMENSION  X(N),SX(30),ISX(30) 
DIMENSION  CV(30) 

THIS  SUBROUTINE  DOES  THE  Rll  TEST  FOR  OUTLIERS  AT  THE  .05  CRITICAL 
TEST  LEVEL. 

THE  CV  VALUES  ARE  THE  CRITICAL  VALUES  FROM  TABLE  V. 

DATA  CV/3*0,. 955,. 807,. 689,. 610,. 554,. 512,. 477,. 450,. 428,. 410, 
1  .395,. 381,. 369,. 359,. 349,. 341,. 334,. 327,. 320,. 314,. 309,. 304, 
1  .299, .295, .291, .287, .283/ 

INPUT  VALUES  ARE:   X  -  THE  DATA  TO  BE  TESTED 

N  -  NUMBER  OF  DATA  POINTS 

OUTPUT  VALUES  ARE:   NOUTB  -  FLAG  FOR  SMALLEST  VALUE 

0  IF  OK,  1  IF  NOT 
IB    -  LOCATION  IN  X  ARRAY 
NOUTT  -  FLAG  FOR  LARGEST  VALUE 
0  IF  OK,  1  IF  NOT 
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C  IT    -  LOCATION  IN  X  ARRAY 

C  IER   -  ERROR  RETURN  INDICATOR 

C  0  IF  ALL  OK 

C  1  IF  N  OUT  OF  RANGE  5  TO  30,  INCLUSIVE 

C 

C 

C     CHECK  FOR  SUITABLE  N 

C 

IER  =  0 

IF(N.LT.4.OR.N.GT.30)THEN 
IER  =  1 
RETURN 
ENDIF 

DO  100  1=1, N 
SX(I)  =  X(I) 
ISX(I)  =  I 
100  CONTINUE 
C 

C     SORT  INTO  INCREASING  ORDER 
C 

CALL  SHSORT(SX,ISX,N) 
RB  =  0. 

IF(SX(N-1)-SX(1)  .NE.  0.)  RB  =  (SX(2)  -  SX( 1 ) )/(SX(N-l )  -  SX(1)) 
NOUTB  =  0 

IF(RB.GT.CV(N))THEN 
NOUTB  =  1 
IB  =  ISX(l) 
ENDIF 
RT  =  0. 

IF(SX(N)-SX(2)  .NE.  0.)  RT  =  (SX(N)  -  SX(N-1 ) )/(SX(N)  -  SX(2)) 
NOUTT  =  0 

IF(RT.GT.CV(N))THEN 
NOUTT  =  1 
IT  =  ISX(N) 
ENDIF 
RETURN 
END 

C    Included  for  completeness  only  -  from  NONIMSL  library  at  NPS 


C  

c 

C  A.  IDENTIFICATION: 

C  TITLE:      NUMERICAL  SORT 

C  ID:         Ml-NPG-SHSORT  (F-IV) 

C  PROGRAMMER:  R.  BRUNELL 

C  DATE:       MARCH  1968 

C  MODIFIED:    DEC.  1973   BY  L.  NOLAN 

C 

C  B.  PURPOSE: 

C  TO  SORT,  IN  ASCENDING  ORDER,  AN  ARRAY  OF  SINGLE  PRECISION  REAL 

C  NUMBERS  BY  THE  METHOD  OF  SHELL,  AND  TO  PRODUCE  AN  ARRAY  OF  INDEXES 

C  SO  USER  CAN  RE-ORDER  OTHER  CORRESPONDING  INFORMATION  ACCORDING  TO 
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C     ASCENDING  VALUES  OF  "A". 

C 

C  C.  USAGE: 

C     1.  CALLING  STATEMENT: 

C        CALL  SHSORT(A,KEY,N) 

C     2.  ARGUMENTS: 

C         A  -  ARRAY  OF  NUMBERS  TO  BE  SORTED.  THIS  ARRAY  IS  SORTED 

C  (RE-ORDERED)  BY  "SHSORT". 

C        KEY  -  ARRAY,  DIMENSIONED  AT  LEAST  N  IN  CALLING  PROGRAM,  TO  BE 

C  FILLED  BY  USER  WITH  INTEGERS  FROM  1  TO  N.   AFTER  EXIT 

C  FROM  SHSORT,  KEY(l)  WILL  CONTAIN  THE  ORIGINAL  INDEX  OF 

C  THE  SMALLEST  ELEMENT  OF  "A";  KEY(2)  WILL  CONTAIN  THE 

C  ORIGINAL  INDEX  OF  THE  NEXT-TO-SMALLEST  ELEMENT  OF  "A"; 

C  ETC.   KEY(N)  WILL  CONTAIN  THE  ORIGINAL  INDEX  OF  THE 

C  LARGEST  ELEMENT  OF  "A". 

C  N  -  NUMBER  OF  MEMBERS  IN  ARRAYS  "A"  AND  "KEY". 

C 

C  D.  REFERENCES: 

C     1.  "ALGORITHM  201,  SHELLSORT",  BOOTHROYD,  J.,  "COMMUNICATIONS  OF 

C         ACM",  VOL  6,  NO  8,  AUGUST  1963,  P. 445. 

C     2.  "CERTIFICATION  OF  ALGORITHM  201",  BATTY, M. A.,  "COMMUNICATIONS 

C         OF  ACM",  VOL  7,  NO  6 ,  JUNE  1964,  P. 349. 

C 

SUBROUTINE  SHSORT ( A, KEY, N) 

DIMENSION  A(N),KEY(N) 

Ml  =  l 
6  M1=M1*2 

IF  (Ml  .LE.  N)   GO  TO  6 

Ml=Ml/2-l 

MM=MAX0(M1/2,1) 

GO  TO  21 

20  MM=MM/2 

IF  (MM  .LE.  0)   GO  TO  100 

21  K=N-MM 

22  DO  1  J=1,K 
II=J 

11  IM=II+MM 

IF  (A(IM)  .GE.  A( II  ) )   GO  TO  1 
TEMP=A(II) 
IT=KEY(II) 
A(II)=A(IM) 
KEY(II)=KEY(IM) 
A(IM)=TEMP 
KEY(IM)=IT 
II=II-MM 

IF  (II  .GT.  0)   GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 
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